/**
 SofferBound.cxx

 \file
 Implementation of class SofferBound.

 \author Thomas Burton 
 \date 3/13/12
 \copyright 2012 BNL. All rights reserved.
*/

#include "tppmc/SofferBound.h"

#include <cmath>

#include <boost/foreach.hpp>

#ifdef USING_LHAPDF
#include <LHAPDF/LHAPDF.h>
#else
#include "tppmc/Enumerations.h" // Contains copy of LHAPDF::Flavour enum
#endif

#include "tppmc/Dssv.h" // For default helicity distribution
#include "tppmc/Grv98Pdf.h" // For default unpolarised distribution without LHAPDF
#include "tppmc/LhapdfDistribution.h" // For default LHAPDF unpolarised distribution
#include "tppmc/parameters.h"
#include "tppmc/Particle.h"

namespace tppmc {
   namespace transversity {

      SofferBound::SofferBound()
      : mEnableSea(false) {
      }

      SofferBound::SofferBound(boost::shared_ptr<PartonDistribution> f,
                               boost::shared_ptr<PartonDistribution> g,
                               bool enableSea)
      : mUnpolarised(f)
      , mHelicity(g)
      , mEnableSea(enableSea) {
      }

      SofferBound::~SofferBound() {
      }

      int SofferBound::Initialise() {
         // Initialise default PDFs if they were not set in the constructor
         typedef boost::shared_ptr<PartonDistribution> PartPtr;
         if(not mUnpolarised.get()) {
#ifdef USING_LHAPDF
            mUnpolarised = PartPtr(new LhapdfDistribution(20250, 0));
#else
            mUnpolarised = PartPtr(new Grv98Pdf(Grv98Pdf::LO));
#endif
         } // if
         if(not mHelicity.get()) {
            mHelicity = PartPtr(new Dssv());
         } // if
         mUnpolarised->Initialise();
         mHelicity->Initialise();
         mPartons = Lhapdf(); // Requires access to PDFs
         return 0;
      }
      
      void SofferBound::Evaluate(const Parton& parton, double Q2,
                                 PartonContent& values) const {
         // Evaluate the unpolarised, f, and helicity, g, distributions.
         // If |g| > f, force |g| = f to obey positivity.
         PartonContent f(0.), g(0.);
         mUnpolarised->Evaluate(parton, Q2, f);
         mHelicity->Evaluate(parton, Q2, g);
         // Set the helicity value of each quark we are supposed to compute
         values = PartonContent(0.); // Initialise output to zeros
         BOOST_FOREACH(int p, mPartons) {
            // Ensure helicity observes positivity for this unpolarised PDF
            if(fabs(g[p]) > f[p]) {
               if(g[p] > f[p]) g[p] = f[p];
               else if(g[p] < f[p]) g[p] = -f[p];
            } // if
            values[p] = 0.5 * (f[p] + g[p]);
         } // BOOST_FOREACH
         // d quark transversity is negative
         values[LHAPDF::DOWN] = -values[LHAPDF::DOWN];
      }

      bool SofferBound::Supports(const Parton& parton) const {
         return mUnpolarised->Supports(parton) and mHelicity->Supports(parton);
      }

      std::string SofferBound::Name() const {
         return "Soffer bound: " + mUnpolarised->Name() + " & " +
         mHelicity->Name();
      }

      std::set<int> SofferBound::Lhapdf() const {
         // Assemble the potential list of partons to check, by PDG code
         std::set<int> pdg;
         pdg.insert(2); // up
         pdg.insert(1); // down
         if(mEnableSea) {
            int sea[8] = { -5, -4, -3, -2, -1, 3, 4, 5 };
            BOOST_FOREACH(int q, sea) {
               pdg.insert(q);
            } // BOOST_FOREACH
         } // if
         // Select the potentials that are supported by both the unpolarised
         // and helicity distributions
         std::set<int> types;
         BOOST_FOREACH(int q, pdg) {
            Parton parton;
            parton.type = q;
            if(Supports(parton)) {
               types.insert(parton.Lhapdf());
            } // if
         } // BOOST_FOREACH
         return types;
      }

      Int_t	SofferBound::Write(const char*, Int_t, Int_t) {
         return 0;
      }

      Int_t	SofferBound::Write(const char*, Int_t, Int_t) const {
         return 0;
      }
   } // namespace transversity
} // namespace tppmc
